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Numerical  Qua.drature  over  a,  Rectangular  Domain 
in  Two  or  More  Dimensions 

By  J.C.P.  Miller 
III  Quadrature  of  a.  Harmonic  Integrand 

Introductory 

1.  In  Note  1(1),  §5,    formula  (B),<^7,  formula  (B'),  and  in  §  9,  also  in  Note  Il(l) 
in  several  places,  we  have  seen  how  the  error  term  is  very  much  reduced  if  the 
integrand  f(x,y)  is  a  harmonic  function,  that  is,  if w  f  -   0.   In  this  note  we 
pursue  further  this  special  case,  in  which  especially  high  accuracy  is  attainable 
with  few  points. 

It  may  not  be  often  that  the  integrand  will  have  this  special  form,  but  it 
seems  worthwhile  to  develop  a  few  of  the  interesting  formulae.  We  start  by 
obtaining  expansions  for  n  variables,  a.nd  more  extensive  ones  for  two  variables, 
and  then  obta.in  and  consider  special  quadrature  formulae. 

2.  Expansions.  As  in  11(1),  §2,  ve  develop  f(^,  V...,xJ  as  a.  Taylor  series 
in  even  powers  of  ea.ch  of  the  variables  Xj,.   Then,  using  v  f n  =  0  whenever  it 

is  applicable,  we  obtain 


rj  =   l/(2h)n  =   (2h)~n/j^  J      f(Xl,x2,...,Xn)    dx1, 


dx_ . . .  dx 
2  n 


.  f    +^    itj3f    +*L    i6^6      +h!(l6    8      192^ 
o  +  5.*     3^0?!       3  J    0      9.'  v  5  **         5   <*-  '  o 

h10  ,128  _>      6      1280      10. 


(2.1^  +ttt(—  rf  a  +±*t<P  )fo 

2r64      12      U096      12      61181^4  ^8      707384     /L2. 


k12 
+  13. 


V 

where,    as  before,    extended, 

6  8 

(2:2)  nr\  -  Z  jV,,        Cr\  -  s    /f°     a  Jl\  -  z      g    g/°2     g 

»Xr3xs  ^r^xt  8xr  8xs  3\  *xu 

etc.,   the  summations  extending  over  all  possible   combinations  of  r,    s,    t, ...with 
no  two   equal. 
(l)  J.C.P.Miller,M.T.A.C,    in  the  press. 
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We  have  likewise,  labelling  the  symmetrical  sets  of  points  as  in  Note  II, 
the  expansions  for  sums  of  values  of  f  over  the  sets 


(2.31)  0  fQ 


(2.32)     a(a)  2nfQ     -  ***-  tftQ  +  *£-  rftQ  +  ^(^^ 

-10h10a10     '■  -1212 


10 


B»V/V0  -  ^W^c^-e^/V--- 


(2.33)   P   (b)   2n(n-l)f0-^Mn-10K?f0    -  i2^  U-i6)  J6f0 


IP      IP 

>  ^,   f2(n-34)^12-3(n+362)C7L2-6(n-728)^<^-f  6.(n-102^12J?  f Q+   . . . 
(2.3^)aT(c,d)  Mn-l)fQ  -  ^-{(n-l)(Ad4)-6c2d2j^f0+  i|^/(n-l)(c6+d6)-15c2d2(c2+d2)))^fD 


8h8rr/     lV  8   .8v  Oo  2^2,  k    Ik  „  kMJS 
+  757—  U(n-l)(c  +d  )-28c  d  (c  +d  )+70c  d  jtr 


-2[(n-l)(c8+d8)-28c2a2(cVdU)-70oV^]   fQ 

+     ..  . 


(B.35)<C(e)  |n(n-l)(n-2)f0-  ^|-(n-2)(n-7)^f0  +  ^2|r2-  (n2-33n  +  122)C76f0 


Q      O 

^{(n-2)(n+13)T38   -2(n2-129n  +  10^^  f( 


We  recall  that  0  is  the  origin,  or  centre  of  the  square,  a  (a)  includes  all  points 
with  one  coordinate  +  ah  and  the  rest  zero,  P(b)  has  two  coordinates  each 
independently  +  bh  and  the  rest  zero,  V  ( c, d)  has  one  coordinate  +  ch,  another  +  dh 
and  the  rest  zero  a.nd  finally£(e)  has  three  coordinates  each  independently  +  eh 
with  the  rest  zero. 

3.  Expansions  over  a.  square  are  simpler  since.!/  f„,  ^~   f_  etc.  vanish.   They  can 
"be  obtained  by  analysis  with  the  detached  operators  -  in  particular /vj  we  proceed 
to  obtain  expansions  with  general  terms . 
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If  F(z)  =  u  +  iv  is  a  function  of  a  complex  variable  z  =  x  +  iy  then  both 

2i    2i 

u  and  v  are  harmonic  functions  satisfying  D  cj)  +  D  9  =  0.   Likewise,  if  u  is  a 

harmonic  function,  it  can  he  shown  that  v  exists  such  that  u  +  iv  is  a  function 
of  a  complex  variable.  We  then  have 

D  F  =  IF'  =  iD  F 
y         n 


D  D 

x  y 


■**-< 


ID  =  -i  D 
x      y 


and 

(3-D 

In  order  to  develop  expansions  we  therefore  substitute 

(3.2)  dx  =  i'^lr   Dy  =  i1/^ 

Consider,  firstly 

_   /-h  s*l  .       hD   -hD    hD   -hD 

(3.3)  J  =  (2h)*2  /    /   f(x,y)dx  dy  =      X     (e  X-e   x)(e  y-e   y)f 

-h  u-h  kh     D  D  ^ 

x  y 

The  operator  is 

sinh  h®x  3inh  hPy  sinh  i     '  ~h^"sinh  i  '      h  %T 


h     D       D 

x      y 


(3.10 


h2^2 

_  ,.1/2     .-1/2^,  ^  .    ,.1/2  .-1/2^^ 

1     cosh  (1  /    +   1     /    )h73   -   cosh  (1  '    -1     '    )hFT 

=  2  h£S2 

1     cosh  y2  h     -cos   >|2  h'fcr 

"2  , 2  _2 

h     TJ 

whence 

(1   *\        t        I*  2V4    ^  25h8     -JB  22r+Vr         Lr  ,    . 


ahD  -ahD  a,bD  -ahD 

Likewise/' Z  f(x   ,    y   )   =   (e  +  e         X+  e       y+  e 

=  2   (cosh  ah  D     +   cosh  ah  D   )    f_ 
v  x  y       0 


)    tr 


=  2   (cosh 


i"1/2  ahJfr+   cosh  l1'2  ah*5)   tr 


(3.6) 


=  h   cosh 


-  Mi  - 


ah 


^<?D    ?r 


2"^       0 


a     h     s&  m-       a, 


$ 


8U8      q  ^rv>r    ). 


]f, 


-  k  - 


and 


Z       f(x_,y_)   =  k  cosh  bh  D     cosh  bh  D     f_ 
fp(b)  P     P  X  y     ° 


=  4  cosh  i-^bhJbcosh  W2  bh  Sb  f 

=  2  (cosh  bhN/2~£T+  u       bh  \/~2^0    f 

,  r-,      2     b     h  <=>.4       2     b     h    ^.d 

=  k  [1+  n T3     +  5! W      +  . 


IT] 


8: 


-2r.kr,kr       \ 

b     h       ^+..]f 


V^T 


o 


We  shall  not  use  all  the  expansions  given  above  in  the  present  note,  but 
it  seems  useful  to  set  out  the  collected  results  for  future  use, 

k.      Lattice-point  formulae  over  a,  Square.   We  consider  first  formulae  in  two 
variables,  and  start  with  9  points,  putting  a  =  b  =  1  and  using  the  sets  0,  Qi(l), 
3(l).  We  write 


(4.1)  J  =  l/4h  =  AQ  f(0,0)  +  Z  Aq  f  (Xq,  ya)   +  Z  Ap  f  (Xp,yp) 
using  (x  y  )  etc.  as  typical  sets  of  coordinates. 

Using  (3-5)  -  (3-7),  we  equal  coefficients  of  Tj   fn,  r  =  0(l)2. 


This  gives 


(*.2) 


AA  +  U  +  U  =  1 
0     a     3 

k 

-  h  A     +  16  A  =  - 

a     3   15 
+  4  A  +  6,  AQ  =  i| 


12 


with  correction  term  C  =  -  ( -k  A     +  256  Aa  -  -=-=-)  — 

OJ  p         yi         -L2 


r^Sfo 


We  obtain  the   formula 


(fc.3) 

^ 

-32 

7 

-32 

1000 

-32 

7 

-32 

7 

■7  900 


14.1.      4  4--  19^2  h^  «_L2   - 

with  main  correction  term-  ^>  ^^T  ~    ^ 


This  formula  is  remarkably  good.  With  the  example  of  Part  I,  we  have, 

2 
writing  J '  -  h  J 

1    /1.2    /d.2  x 

J'  =  r-    /  sin  x  sinh  y  dx  dy  =  t-  (l-cos  1.2)(cosh  1.2  -l) 


Formula,  (4.3)  gives 


,12 


=  0.12922  70590  73675  11602 

JT  =  0.12922  70590  72834  11029 
12/ 


withE=  -0.-0  84l  00573  and  C  =  +0.0  8^1  01633. 
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5.   Five-point  formulae.   The  high  precision  of  (^.3)  suggests  that  formula.e 
of  lessor  precision,  with  fewer  points,  may  be  useful.  We  use  the  first  two 
of  (^.2)  and  take  one  of  A„,  A  ,  A  to  be  zero|(i)  A0  =  0  gives  an  eight -point 
formula  with  relatively  poor  precision. 


(5.1) 

19 

56 

19 

56 

0 

56 

19 

56 

19 

(ii)   A 

=  0 

gives 

(5.2) 

l 

0 

l 

0 

56 

0 

l 

0 

i 

300 


,8     g 
with  main  correction  term  -k-0  frrTr    f 


0 


6o 


with  main  correction  term  - 


32  h8  _8 
-    9T   ^    fo 


(iii)  A.  =  0  gives 
P 


(5.3) 


0 

-1 

0 

-1 

19 

-l 

0 

-l 

0 

7  15 


with  main  correction  term  +  —  — 


5  9 


r*8 


0 


7     8 

We  observe  tha.t  (5.2)  and  (5«3)  combine  in  the  proportions  ir=   *  r-=-  to  give  (^.3).? 

7  k  1J   J-5 

though  without  an  error  estimated   Likewise  *■  X  (5«2)  -  -~   >(  (5-3)  gives  (B)  of 

lip-5  1.8   a 
Note  I,  and  an  estimate  for  the  correction,  namely  -  ±±£   £_  *^)°  ^     when  (x  v)  is 

harmonic . 

Another  combination,  that  of  (5-2)  and  (5.3)  in  equal  proportions  gives  a 
small  correction  term: 


(5.M 

l 

-k 

1 

-k 

132 

-k 

l 

-k 

L 

Again  k 

K   (5.2) 

-   3 

*(5.3 

)  § 

(5.5) 

1 

3 

l 

3 

-1 

3 

l 

3 

l 

*  120 


with  main  correction  term-  p  ttt  <:W  f~ 

5  9»      0 


r  15 


with  main  correction  term  - 


212 


8 


Sr*8 
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Evidently  (4.3)  is  most  precise,  "but  simultaneous  use  of  (5.2)  and  (5.3) 
gives  an  idea,  of  the  precision  attained  and  readily  yield  the  better  result  if 
desired.   Formula  (5«5)  might  be  helpful  with  desk  computing  but  (5.1)  has 
little  to  recommend  it. 

Numerical  results  for  some  of  the  formulae  using  the  example  of  £  4  are  as 
follows : 


Formula. 

Result  J' 

1010/ 

E 

1010  K  c 

(5.1) 

0.12922  72986 

+2395 

-2396 

(5.2) 

0.12922  70974 

+  383 

-  383 

(5.3) 

0.12922  70255 

-  336 

+  335 

I  (B) 

0.12922  71932 

+1341 

+1342 

(5  JO 

0.12922  70615 

+  24 

-  2k 

General  n; 

2 
2n  +  1  points. 

We 

consider 

now 

the 

n-dimensio 

n  >  3,  using  lattice  points  0,  Q:(l),  |3(l).   In  this  case  the  term  in  J  fn  is 

o  0 

relevant,  and  the  4L  f n   term  will  appear  in  the  error,  except  when  n  =  3. 

k        *J& 
We  equate  coefficients  of  f~,  sT   f_, u   f  in  the  expansions  resulting 

from  use  of  (2.1),  (2.3l)-(2.33)  in  (4.1).  We  obtain 

(6.1)      ( AQ   +  2n  AQ  +  2n(n-l)  Ap  =  1 

-  k  A     -   8(n-4)  AQ  =  ~ 
a    v  '     3   15 

+  6  Aa  +  I2(n-l6)  Ap  =  ~ 

o  o 

while  0  .   -{k  Aa+  8(n+6)Ap-  §§r*\  +K  ^(n-6k)A^+  iff  }   jjj.  s}   f Q 

These  yield 

((,  o\  A    -6ln2  +  931n  +  3780        6ln  -496  -6l 

K°'d)  A0  "      3780  a  "   3780       Ap  *   T55o 

with  „   II98   h8  ^8  _    3619   h8  <±Q  _ 

c  -  -95^  BT*7  f 0  +  -315   BTXfo 

In  particular 

^  a^        a     a    12048      .       626       .        61 
(6'33)    n  =  3     A0  =  "75^      Aa  =  "  7560       AP  =  "  75SO 


(6  34)    n  -  4     A  -  i32S§      A St      A  =  -  — & 

^o.j4;    n  -  4     aq   -^q      Aa     ^q       Ap     y^Q 


Aa  "   7560 

V- 

A         5°U 

a  "   756o 

V- 

-  7  - 
tti   ^^        c     «    1382°       a      382      .      6l 

(6  36)     n-6     A  =^£        A  _._260       A-__4i 
lD,jDj    n  _  D    A0   T5S0      Aa    75^0      V  75^0 

As  a  numerical  illustration  for  n  =  3  we  consider 

1    1    r  1  r1  C  1  3  5 

J  =  7T  I  =  7T  cos  jf  x  cos  y  cosh  f  z  dx  dy  dz 

J  .1  J-i  Jjl 

=   ~  sin  r  sin  1  sinh  r  =  O.98OO827 
15      4  4  ' 

Formula  (6.33)  gives  J  =  0*9799734  WithE=  -0.0000109  and  C  =  +0.0000110. 

This  result  is  less  spectacular  than  that  of  Sk,  for  these  reasons: 

i)  In  &4,  h  =  0.6.,  here  h  =  1,  a.nd  the  correction  term  in  (4.3)  contains  a  high 

power  of  h. 

8  12 

ii)  The  correction  term  in  (6.2)  is  of  order  h  ,  that  in  (4.3)  is  of  order  h 

iii)  The  higher  the  number  of  dimensions.,  the  more  individual  terms  there  are 

in  cD  f }  JD^  f<>  etc.   In  (4.3)  there  is  only  one  term  in<£)  f3    in  (6.33) 

o 

there  are  9  in^  f* 

iv)  The  effect  of  larger  interval  h  is  enhanced  by  the  use  of  the  factor 

5  5 

x-y    which  exceeds  unity.,  in  cosh  t-  Z;  this  is  only  partially  balanced  by 

3 

the  factor  cos  j~   x. 

In  spite  of  these  points ,   the  formula  (6.2)  seems  a  good  one. 

7.   Quadrature  over  a.  square^,  specially  chosen  points.   Since  the  expansions 
of  §  3  involve   only  cross-differences  2)  f_,  it  appears  likely  that  use  of 
sets  of  diagonal  points  f3  will  be  mo."-  profitable  than  attempts  to  use  sets  a. 
It  turns  out  that  sets  0,  a(a),  p(b)  and  0,a(a.)  both  fail  to  give  real  values 
of  a  if  maximum  precision  is  sought.   On  the  other  hand,,  we  can  get  several 
formulae  making  use  of  any  number  of  sets  ^(b^) , -\>  =   0(l)  r,  both  with  and 
without  the  point  0. 

We  start  first  with  r  sets  p  ^bp).,  without  the  point  0.  We  have  to  find 
the  2r  constants  A  3    bj,  satisfying  the  equations 
r 

(T-1)      1   ,  k\     b^(S"l)  =  T2s-l)Os-3)  "  \-V     s  "  1(1)  2r 

obtained  by  substitution  of  (3.5)  —  (3. 7)  in 
(7.2)         J  =  I  Aft  f  (+  b.  h,  +  b.  h) 


-  8  - 


and  equating  the  coefficients  of  the  first  2r  coefficients  of  - 

powers  of  k   have  been  cancelled. 

k 
By  familiar  arguments,  the  bi,  are  rqots  of  the  equation 

2 


Sundry 


(7.3) 


x 

c 


c    c  -  c  _ 

r     r-1   r-2 


r+1 


'2r-l 


=  0 


These  are  the  orthogonal  polynomials  corresponding  to  the  weight  function 

1   -3/^-   -l/2 
w(x)  =  ;r  (x  '   -x  '  )  for  range  0  <  x  <  1.   The  first  two  are 


(7.10 


I5x  -1=0 

8l9x2  -  k3Qx  +   11  =  0 


(7.5) 


The  main  correction  term  is  obtained  from  the  next  power  of  S)  and  yields 

8r  x  M,  8r  <s8r 


C  ..  (C_  -  Z     k  JL  b°r  )   2^ 


2r 


If  the  point  0  is  included,  our  equations  (7»l)  are  replaced  by 


A„  +  klA„     =  1 


0 


3% 


(7.6) 


\>=i' 


£     4  A„     b 


biL       %     *     =   (2s+l)(Wl)     =   S 


=   C    ,      s  =  1(1)   2r 


k 


and  the  b*>  are  roots  of  the  equation 

2 


1 

X 

X 

cl 

C2 

C3 

°2 

C3 

ck 

°r+l 

°r+2 

°r+3 

'r+1 


'r+2 


'2r 


=  0 
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which  are  the  orthogonal  polynomials  for  the  weight  function  w(x)=  ^(x  '  -x  ' 
for  the  range  0  <  x  <  1.  The  first  two  are 

3x  -  1  =  0 
(7-8)  2 

I7017x  -  1365OX  +  1745  =  0 

The  main  correction  term  is  this  turns. 

r         o  .   02r+l  ,4r+2   i,  „ 

(T.9)       c  .  (c2r+1  -£  4  ^  ^)  ^f^  »^2 

In  each  case  the  coefficients  Ac  may  be  computed  hy  standard  methods. 

p 

Kr 

8.  Formulae  for  r  =  1.  These  ha,ve  4  and  5  points  respectively 

8^8 


(8.1)  Ap  =  ¥   b  =  15  C  =  225   "ST"* 


(8.2)  An  = 


4 


o 

12  iv  12 


a    1      h  q-ia    p  -  2816  h^  ^  f 

0  "  5     P  "  20      °  "  5  U  ~  122H5      121   J° 

Written  out  in  full: 

(8.3)  j -J  1  -\ (fdr^ir^j^-is-^ir^Hfds-^-is-^M-is^s-^j 

(8.4)  j.  I  1=  |f(o,0)^|f(3-1A3-^),f(-3-^,3-^)+f(3-lA,-3-^).f(-3-^-31/i)} 


Asa  numerical  test  use 


t   i  t   1  r x  r 


cos  x  cosh  y  dy  =  sinl  sinhl 
=  O.98889  77057  62865 


Formula  (8.3)  gives  0. 98889  06525 

with      E  =  0.00000  70533  and  C  =  +0.00000  705+7 
o^tuL      formula  (8.4)  gives  0. 98889  77062  +1358 

with  E,.=  +0„094  78^93  and  C  =  -0. 0^785^3. 

9°  Formulae  for  r  =  2.   These  have  8  and  9  points  respectively 

(9.1)  b_  =  0.40316  26030  59346  89754     AQ  =  0.22912  30654  28169  97222 

b2  =  0.84439  75319  23478  74713    Ap  =  0.02087  69345  71830  02778 

54592     256  ,16  0.16 
with  mam  correction  tem—r^nr  X  TZT     n  *D        f 
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(9.2)     bQ  =  0  AQ  =  0.69521  80834  12925  81989 

bn  =  0.63205  02078  I8796  99524  Ac  =  0.06686  42185  46105  38-162 

1  f3]_ 

b2  =  0.89531  6379I  24106  97730  A^=  0.00993  12606  00663  16340 

,4.x.                            +•           16832  v  1024  ,20  q,20  . 

with  main  correction  term   w   qq  a  pn,  h   JO    f 


For  J  =  i 


r  1  ri 

cos  x  cosh  y  d/  dy 

-1    v/-l 


formula  (9.1)  gives  0. 98889  77057  62853  38396 

withE=  -0.013  II71243  and  C  =  +0.0131171555 

while  formula  (9.2)  gives  0. 98889  77057  62865  09647 

with  E  =  +0.0199  and  C=  0.0^90 

■With  formula  (9«2)  we  find  approximately 
0.82447  37090  77903  16756  for 

1  r  2  r 2 

-r-?-  I     cos  x  cosh  y  dx  dy  =  sin  2  sinh  2 

=  0.82447  37090  77809  15433 
with  E  =  0.0139401323  and  C  =  0.01394o6250. 

These  formulae  clearly  have  high  precision,  even  with  considerable  values 
of  h. 

10.   Quadrature  over  a.  Cube;  specially  chosen  points.   The  search  for  such 
formula.e  is  more  difficult  in  3  or  more  dimensions.   It  seems  that  one  or  more 
extra  available  constants  are  needed  in  order  to  obtain  real  points.  We  shall 
not  pursue  this,  but  give  one  simple  formula  for  three  dimensions. 

We  find  nothing  convenient  by  use  of  points  a(a),  with  or  without  0}    likewise 
0  with  (3(b)  fails  to  give  real  points.  We  can,  however,  use  12  points  p(b) 
alone.  We  have  then  to  satisfy 

2n(n-l)AQ  =  1 

4        4 
8b  (4-n)A  =  — ,  where  n  =  3- 

with  main  correction  term. 

l0  H  21'  6.'  ~      "  —~'~-  "  ^   *o 


This  yields   b  =  (2/5)  '   =  0,79527  07287  67051   A  =  l/l2 

c  =  (156  b6  a  +  |f)  |j. a6  f0  -  0.00562  h6  a6  fc 


-  11  - 

With  the  example  of  £6,  with  integrand  cos  r-  x  cos  y  cosh  z 
(10. 1)  gives   J  =  0.97519  with  £  =   -0.00489  and  C  =  +0.00^. 

-.6 

The  only  formula  found  that  allows  for  the  term  sj      f  and  ha.s  a.n  error 
o  0 

of  order  h  is  (6.33).?  which  needs 9  points.   It  is  evident  that  further  search 
is  needed. 
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